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Abstract 



^ . A high-order scheme for direct numerical simulations of turbulent combustion is 

I ! discussed. Its implementation in the massively parallel and publicly available Pen- 

2 ! ciL Code is validated with the focus on hydrogen combustion. Ignition delay times 

'^ ' (OD) and laminar flame velocities (ID) are calculated and compared with results 

, ^ , • from the commercially available Chemkin code. The scheme is verified to be fifth 

order in space. Upon doubling the resolution, a 32-fold increase in the accuracy 
of the flame front is demonstrated. Finally, also turbulent and spherical flame 
front velocities are calculated and the implementation of the non-reflecting so- 
called Navier-Stokes Characteristic Boundary Condition is validated in all three 
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i/^ • directions. 
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1. Introduction 

^ . Modeling of turbulence is one of the largest research areas within flow me- 

chanics. Turbulent combustion inherits all the properties of non-reacting turbu- 
lent flow. The most important addition is linked to the highly nonlinear reaction 
processes, and models for this are called combustion models. Two additional 
challenges in turbulent combustion are the very sharp changes in density and dif- 
ferential diff"usion of mass and heat. 

For combustion processes it is crucial to be able to simulate the mixing of 
the combustible species correctly. Traditionally this has been done by means of 
mixing models in Reynolds Averaged Navier Stokes (RANS) codes by combining, 
e.g., the k-e turbulence model and the eddy dissipation concept (or EDC) "mixing" 
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model 1231 |2(a 12211 . or in Large Eddy Simulation (LES) 12711 where a sub-grid 
model is used both for the turbulence and for the scalar mixing. There are however 
major and still unresolved problems related to modelling of what happens on the 
very smallest scales with these methods. 

Several RANS codes with detailed chemistry are commercially available [jg, 
I12I1 . and there are a huge number of these codes found as in-house codes at dif- 
ferent academic institutions and in many industrial departments around the world. 
There are also freely available open- source RANS codes with detailed chemistry 



[|19ll . The reason for the popularity of RANS is its low demand on computational 
resources. Because of this RANS has, for decades, been the most used type of 
code for industrial purposes. 

Nevertheless, also LES has increased in popularity during the last years, and 
this has led, for example, to the inclusion of a LES module in [12]. Most LES 
codes for combustion today are, however, in-house codes owned by different aca- 
demic institutions. 

The most accurate way of simulating turbulent combustion is to use Direct 
Numerical Simulation (DNS) |127D instead of RANS or LES. In DNS one resolves 
the full range of time and length scales of both the turbulence and the chemistry 
(using accurate high-order numerical methods for computational efficiency). The 
problem with DNS is however that it is very resource demanding, both on CPU- 
hours and memory. 

In this paper we present the implementation of a detailed chemistry module 
in a finite-difference code flj for compressible hydrodynamic flows. The code 
advances the equations in non-conservative form. The degree of conservation of 
mass, momentum and energy can then be used to assess the accuracy of the solu- 
tion. The code uses six-order centered finite differences. For turbulence calcula- 
tion we normally use the RK3-2N scheme of [|32ll for the time advancement [5j. 



This scheme is of Runge-Kutta type, third order, and it uses only two chunks of 
memory for each dependent variable. For hydrodynamic calculations, the lengths 
of the time step is calculated based on a number of constraints involving maxi- 
mum values of velocity, viscosity, and other quantities on the right-hand sides of 
the evolution equations. In some cases we use instead a fifth-oder Runge-Kutta- 
Fehlberg scheme with an automatic adaptive time step, subject to the aforemen- 
tioned hydrodynamic constraints. However, in many cases we found it advanta- 
geous to use a fixed time step whose length is estimated based on earlier trial runs 
with an automatically calculated time step. 

On a typical processor, the cache memory between the CPU and the RAM is 
not big enough to hold full three-dimensional data arrays. Therefore, the Pencil 



Code has been designed to evaluate first all the terms on the right-hand sides of the 
evolution equations along a one-dimensional subset (pencil) before going to the 
next pencil. This implies that all derived quantities exist only along pencils. Only 
in exceptional cases do we allocate full three-dimensional arrays to keep derived 
quantities in memory. However, most of the time, multiple operations including 
the calculation of derivatives is performed without using intermediate storage. 

As far as we are aware, no open source high-order DNS code with detailed 
chemistry is currently available. The amount of man-hours for implementing a 
fully parallelized DNS code with detailed chemistry is enormous. It is therefore 
now timely to make such a code available in the public domain and to encourage 
further development by a wider range of scientists. Here we describe the imple- 
mentation of such a scheme in the Pencil Code, which is currently maintained un- 
der the Google Code subversion repository, http : //pencil- code . googlecode .com/ 
The code is highly modular and comes with a large selection of physics modules. 
It is portable to all commonly used architectures using Unix or Linux operating 
systems. The code is well documented and independent of external libraries and 
any third party licenses. All parts of the code, including the current chemistry im- 
plementation, is therefore explicitly open source code. In particular, there are no 
pre-compiled binary files. Consequently there are no licenses required for running 
any part of the code. It is therefore straightforward to download the full source 
code from the original subversion repository on google-code. The Message Pass- 
ing Interface libraries are needed when running on multiple processors, but all 
parts of the code can also run on a single processor without these libraries. The 
integrity of the code is monitored through the automatic execution of a selection 
of test cases on various platforms at different sites. The detailed history of the 
code with about 14,000 revisions is accessible. 

It should be emphasized that the use of high-order discretization is critical for 
optimizing the accuracy at a given resolution. Doubling the resolution of a 3D 
explicit code require 16 times more CPU time, but this increases the accuracy by 
a factor of 32. In fact, switching to a derivative module with a tenth order scheme 
is straightforward and not significantly more expensive. 

2. The equations 

In this section we present the governing equations together with the required 
constituent relations such as the equation of state and expressions for viscosity, 
diffusivity and conductivity. 



2.1. Governing equations 

The continuity equation is solved in the fonii 

^— ■ 

where D/D? = didt + f/ • V is the advective derivative, p is the density, and U is 
the velocity. The momentum equation is written in the form 

P^ = i(-Vp + F,,)+/, (2) 

ut p 

where p is pressure, / is a volume force (e.g. gravity or a random forcing func- 
tion), 

Fvs = V • (2pyS) (3) 

is the viscous force, where S,y = ^{dUi/dXj + dUj/dXi) - ^S^V ■ U is the trace-less 
rate of strain tensor. The equation for the mass fractions of each species is 

P^ = -^-Jk + cOk, (4) 

where Y is the mass fraction, / is the diffusive flux, co is the reaction rate and 
subscript k refers to species number k. Finally, the energy equation is 

P ml Dt ^ Dt \mt TJ m T pT 

where T is the temperature, Cp is the heat capacity at constant pressure, R is the 
universal gas constant, h is the enthalpy, m is the molar mass, and q is the heat 
flux. The reason for solving for the temperature directly, instead of, e.g., the total 
energy, is to avoid having to find the temperature from the total energy afterwards. 
In this work we use the ideal gas equation state given by 

P = • (6) 

m 

In the following we discuss the detailed expressions for viscosity, reaction rate, 
species difl'usion, thermal conduction, enthalpy and heat capacity. 



2.2. Viscosity 



The viscosity v is the viscosity of the mixture given by OOll 
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where A^^ is the number of species, Vk is the single component viscosity, Xt 
Ykm/nik is the mole fraction of species k, and 
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The single component viscosity is given as [SO 
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where crt is the Lennard- Jones collision diameter,^ is the Boltzmann constant. 
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and f^j^ is the collision integral that is given by [|25|1 
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where ^\^i]\ is the Lennard- Jones collision integral and 
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are the reduced dipole moment and temperature, respectively. In the above equa- 
tions, ek is the Lennard- Jones potential well depth and jj-u is the dipole moment. 

The values of eu, Hk and CTk must be given as input ll24ll . while the Lennard- 
Jones collision integral is represented by 
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where the coefficients a] are found from Table [B 
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Table 1: The a, coefficients are used in Eqs. (fT2] i and ( l20l i and are taken from the paper of lliol] . 



2.3. Reaction rate 

The reaction rate of species k is given by 
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where A^^ is the number of chemical reactions, nik is the molar mass of species 
k, pk is the partial pressure of species k, nk = Pk/in-k is the molar concentration 
of species k, and pk is the density of species k. Furthermore, v[^ and v['^ are the 
stoichiometric coefficients of species k of reaction s on the reactant and product 
side, respectively. The rates of reaction s are given by the Arrhenius expression 



h = BnT"" exp(-Ean/RT\ 



(14) 



where 5„ is the pre-exponential factor, a„ is the temperature exponent, and E^n 
is the activation energy and they are all empirical coefficients that are given by 
the kinetic mechanism. For hydrogen-air combustion, an example of a kinetic 
mechanism is found in rtl7|]. 

2.4. Species diffusion 



The diffusion flux is Jk = pYkVk- Following Bill , the diff"usion velocity, Vk, is 
found by solving 



'^Xp = J]^{y^-Vp) + {Y,-X,)^ + ^J]Y,Ykif,-fk), 



k=i 
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where the Soret effect is neglected. The first term on the right hand side cor- 
responds to ordinary diffusion, the second term is the so called baro-diffusion, 
while the last term is due to unequal body-forces per unit mass among the species. 
Unfortunately the CPU cost of solving Eq. (fT5l) numerically scales as A^^^ for each 
grid point and time-step, and simplifications are therefore required in order to be 
able to run reasonably sized simulations. 

In the mixture averaged approximation the diffusion velocity is expressed as 
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where the body force term has been neglected, D/, is the diffusion coefficient for 
species k 

Dk = ^~^' , (17) 

and Dkj is the binary diffusion coefficient that is given by [il61 

3 ^InklT^/nijk 
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where cTjk = (ctj + cr^)/2 is the reduced collision diameter, mjk is the reduced 
molecular mass for the (j, k) species pair 



mjk = 



tfij + ruk 



(19) 



n*^^'^^* is the collision integral that is given by [llOll 
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where the coefficients a. are also found from Table [B The reduced dipole mo 
ment and the reduced temperature are given by 

- kl and r; = ^, 

2 ejk 



'Jk 



(21) 



respectively, where /j.*^ = fi*iul is the nondimensional 2-species dipole moment, 

ejk = yfejek is the 2-species Lennard- Jones potential, and //^ = Hkl J^kcrl is the 
nondimensional dipole moment. 



2.5. Thermal conduction 
The heat flux is given by 



q = Yj ^kJk - ^Vr, (22) 



where A is the thermal conductivity, which is found from the thermal conductivi- 
ties of the individual species as 
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Here, the individual species conductivities are composed of transitional, rotational 



and vibrational contributions and are given by [12911 
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2.6. Enthalpy and heat capacity 

The enthalpy of the ideal gas mixture can be expressed in terms of isobaric 
specific heat Cp and temperature as 



hi = h1+ CpjdT, h = Y^Yihi, 

J To i^i 



(25) 



where h^ is the enthalpy of formation of species i at temperature Tq. 
To calculate the heat capacity Cp we use a Taylor expansion. 



R ' 



Cp 



= _y«.r-i, (26) 
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where a, are coefficients found in lll5ll . 



3. Scaling in the Pencil Code 

3. 1 . General remarks 

For direct numerical simulations (DNS) it is crucial to have high accuracy. 
This is due to the fact that we are interested in resolving the smallest scales, and 



consequently we can not allow for these scales to be lost due to low accuracy. Fur- 
thermore, for many situations it is important to know the actual Reynolds number 
of the simulation. The Reynolds number is defined as 

ul 
Re = -, (27) 

V 

where u and / are characteristic velocity and length scale, respectively, and v is 
the viscosity. High accuracy is obtained by the use of high order discretization. In 
the Pencil Code, sixth order discretization is normally used [|4i]. However, for the 
density a fifth order upwinding scheme us used. 

3.2. One-step reaction model, R ^ P. 

In order to verify that the code recovers correct scaling, we use sirnplified 
chemistry and compare against known results. Following Doom et al. [9], we 
consider a one-step laminar premixed flame model. The irreversible reaction can 
be presented as R ^> P, where R is the reactant and P is a product. Using the 
approach of Ferziger and Echekki [lllll . we neglect viscous effects, and take p, 
A, Cp and the x-component of the velocity to be constant. Then the system of 
equations takes the form 

dYp dYp _ \ d'^cp . 

dt * dx he dx^ 

dt dx AT pCpT dx^ ' 

where the reaction rate is defined as 

^ ^ I pUlCpA-W - i)(i - Yp) if r > r, 

I otherwise, 

where jS = (Too - To) I {Too - T^), while Tq and Too are the temperature of the un- 
burned and burned gas, respectively, T^ is the critical temperature, Yp is a mass 
fraction of the product, Le = A/{pDCp) is the Lewis number, D is the mass diffu- 
sion coefficient. Taking Le = 1, one obtains the following analytical solution 



^ _ I ± fj '^j^^yA/uj i± A ^ V7, ,^.. 



1 -y8-'exp(jcM) ifjc<0, 

I - /3'^ exp[{l - /3)x/6] otherwise, 

where f = {T - Tq)I{Too - Tq) and 6 = A/{pUCp) is a characteristic thickness. 
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Figure 1 : One-step laminar premixed flame model. Left panel: temperature as a function of x ob- 
tained numerically (solid curve) and analytically (asterisks). Right panel: error of the calculation 
as a function of the mesh spacing 6x is shown by asterisks, and the expected dependence of error 
(proportional to bx^) is indicated by the solid line. 

In the left hand plot of Fig. [1] we compare the numerical results with the ana- 



lytical solution for Tq 
ergg 



-'Y.-\X 



300 K, Too = 2000 K, T, = 440 K, /3 = 1.09, Cp = 10^ 
D = 2 cm^ s"\ p = 5 X 10""^ g cm"^ and 



10^ erg cm-' K"' s-\ 



Ux = 100 cm s'^. It can be seen that there is good agreement between the numer- 
ical and analytical results. To show the high-order spatial accuracy provided by 
the Pencil Code, we obtain the set of solutions for 33, 65, 129, 257, 513 and 1025 
grid points, and compare them pairwise ("33" with "65", "65" with "129" and so 
on). In every pair we compare only the points which are collocated, that is, we do 
the comparison for all the grid points of the coarser grid against half of the grid 
points of the finer grid. The time step is controlled by the chemistry and is here 
fixed at 6t = 10 ** s. The size of the domain is 3 cm. The maximum absolute value 
of the difference between the corresponding solutions in common points is taken 
as the error. In the right-hand panel of Fig. [T]the error as a function of 6x is shown 
by the symbols. One can see that sixth-order accuracy is obtained (see solid line). 

3.3. One-dimensional premixed flame with the Li mechanism 

In this section we study a one-dimensional problem with detailed chemistry. 
We consider hydrogen-air combustion using the Li mechanism [17]. The fresh 
hydrogen-air mixture enters the domain under stoichiometric conditions {Y^^ = 
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Figure 2: Accuracy of calculation as a function of 6x for the Li mechanism. An error of cal- 
culations as a function of the mesh spacing 5x is shown by asterisks, and the solid line is the 
dependence of the error, which is proportional to 5x^ . 



2.4 %, Fo, = 23 % and Fn, = 74.6 %) at a temperature of r„ = 29'&K and a 
pressure of p = 1 atm. To avoid reflection of acoustic waves at the boundaries 
non-reflecting boundary conditions are required. Here the Navier-Stokes Charac- 
teristic Boundary Conditions (NSCBC) [1261 1 ISll have been used. 

To check the spatial accuracy in the case of the Li mechanism, we perform 
numerical experiment as described in Sec. 13.21 for 65, 129, 257, 513, 1025 and 
2049 grid points. The time step is fixed at 6t = 10"'°s. The error as a function of 
6x is presented in Fig.O where one can see that the fifth-order spatial accuracy is 
achieved. The reason we get only fifth order, and not sixth order, is that we are 
using upwinding for the density, which has the efi'ect of decreasing the order of 
the discretization to fifth order [bh. 



4. Validation of the chemistry implementation in the Pencil Code 

In this section the chemistry module will be verified quantitatively by compar- 
ison with the commercially available simulation tool [7]. In order to minimize the 
effect of the fluid flow, and to focus as much as possible on the chemistry, these 
tests have been restricted to zero and one dimensional tests cases. 
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Figure 3: Dependence of the gas temperature on time, computed with the 6-step mechanism (li;- 
amonds), the 8-step mechanism (asterisks), and the Li mechanism {triangles). The numerical 
results are compared with Chemkin for the 6-step mechanism {three doted-dashed curve), for the 
8-step mechanism {doted-dashed curve), and for the Li mechanism {solid curve). 



4.1. Zero-dimensional test: ignition delay 

First a zero-dimensional ignition delay test for different chemical mechanisms 
is studied and compared with the results obtained with Chemkin for the same 
setup. One assumes hydrogen-air combustion in a closed homogeneous reactor at 
constant volume, and consider the 6-step and 8-step mechanisms of um together 
with the Li mechanism [17]. The initial values are p = \ atm for the pressure, 
0=1 for the equivalence ratio, and T = 1200 K for the temperature. As the 
minimum time step varies greatly with the progress of the combustion process the 
time step is here chosen automatically using the adaptive Runge-Kutta-Fehlberg 
method. The results are presented in Fig. |3l where one can see good agreement 
with the Chemkin results. 



4.2. One-dimensional test: laminar flame speed 

Next, we consider a one-dimensional flame front. The cold premixed gas en- 
ters at one end of the domain at given velocity. Inside the domain there is a flame 
front where the fuel is consumed and the temperature increases to the mixture 
flame temperature. The mechanism of P IT"] is used and the inlet values of temper- 
ature, pressure and mixture compositions are the same as described in Sec. 13.31 
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Figure 4: Flame speed velocity as a function of pressure p. Chemkin results are shown by asterisks. 

The inlet velocity is adjusted such that the flame front becomes stationary inside 
the domain. The flame velocity is thus arranged to be equal to the inlet velocity. 

We find that the flame front should be resolved by at least 10 grid points in 
order to ensure a well resolved flame. For a thickness of the flame front of about 
0.01 cm, and a domain of Ax = 0.1 cm, the optimal grid size is found to be 150 
points. The flame speed as a function of pressure is shown in Fig. H] where the 
current results are found to compare well with those of Chemkin. 

5. Three-dimensional flame front simulations 

5.1. Plane flame front 

In this section we study a 3D representation of the initially flat flame front. 
The settings of the problem is similar to that in Section 1431 i.e. initially the tem- 
perature, density and velocity change in the x direction, and are constant in the y 
and z directions. We use periodic boundary conditions in the y and z directions, 
and in the x direction we use inlet and outlet NSCBC boundary conditions on the 
left and right hand sides, respectively, as was done in [llSll . The pressure is p = 1 
bar, the initial gas temperature is T = 750 K, and the inlet velocity is 30 m s ^ 
The unbumed gas mixture has an equivalence ratio of = 0.8. The size of the 
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Run Mechanism Number of species Transport r \ps/Nt/Ng] 



A 


Li 


13 


B 


Li 


13 


C 


No 


13 


D 


No 


13 


E 


No 






Mix- aver. 


71.3 


Oran 


51.7 


Mix- aver. 


42.1 


Oran 


24.3 


Const. 


1.2 



Table 2: Timings r in microseconds per timestep, A^,, per grid point A^^. 



calculated domain is taken to be 0.5 x 0.25 x 0.25 cm^, and the grid size is (128 x 
64 X 64). 

We study both laminar and turbulent regimes. In the laminar regime we check 
that the obtained flame speed is the same as that in the one-dimensional problem. 
In the turbulent case we set the turbulent inlet flux as follows. First, we consider 
an isothermal box with periodic boundary conditions. Initially the density and 
velocity fields in the box are taken to be constant. We use a forcing function in 
Eq. (O similar to that used in Brandenburg [|3|], 

fix, t) = Re{Nfk(r) expiikit) ■ x] + iip(t)}, (32) 

where k{t) is a time-dependent wavevector with k^ = {\k\) being its average value 
that is the chosen to be 1.5 times the minimal wavenumber that fits into the do- 
main, and (p{t) is a random phase. The prefactor A'^ = foCs{kfC^ol6tyi^ is chosen on 
dimensional grounds, c^q is a reference sound speed, and /o is a nondimensional 
factor that it chosen to regulate the strength of the turbulence. 

The simulation is run until the turbulence is statistically stationary. This box 
of statistically stationary isotropic turbulence is then used as the inlet condition for 
the simulation of the turbulent flame front. The values of a two-dimensional slice 
from the box (perpendicular to the main stream) are used as the instantaneous inlet 
velocity, and the slice is changed as a function of time to represent a real inlet. 

For the test case shown here the turbulent intensity is 7 times larger than the 
laminar flame velocity 5l = 10.2 m s"^ We find that for the mean inlet velocity 
35 L the flame is nearly stationary inside the domain. This indicates that the turbu- 
lent flame velocity in this case is around 35 l- However, it is hard to determine the 
turbulent flame speed precisely, because it is difficult to make the flame perfectly 
stationary inside the domain. This is partly because of the fact that between inlet 
and outlet the turbulence is decaying. Far from the inlet the turbulence is weaker 
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Figure 5: From left to right instantaneous mass fractions of H2 and OH are shown. Unburnt 
turbulent gas in injected on the left. 



than close to the inlet, whereas the turbulent flame speed increases with the tur- 
bulent intensity. As a result, the flame which is already far from the inlet tends to 
move even further downstream and the flame brush becomes broader. 

In Fig. [5] one can see that the H2 fuel (on the left hand side of the domain) is all 
consumed over the flame brush. The thickness of the flame brush is of the order 
of half the box length (2.5 mm) and is slightly smaller than the integral scale of 
the turbulence. The mass fraction of OH is shown in the right hand figure. It is 
clearly seen that OH does not burn out after the flame, but due to the very high 
temperatures downstream of the flame front the mass fraction of OH stays rather 
constant. For HO2 the situation is however rather different and it exists only in the 
neighborhood of the reaction zone of the flame (see the left hand figure of Fig. (6]). 
This indicates that HO2 might be used as an indicator of the reaction zone. In the 
right-hand figure the temperature is shown to increase from 750 K to 1984 K, but 
the maximum value will increase even more downstream of the box due to radical 
reconnection. 

In addition, we find that the turbulence is damped behind the flame front, 
and the burnt gas stream looks much more laminar there (not shown here). This 
happens because the values of temperature and hence also viscosity of the burnt 
gas are much larger than those of the unburned mixture. 



15 





Figure 6: As Fig.|5] but for the instantaneous mass fraction of HO2 on the left and temperature on 
the right. 



5.1.1. Timings 

As DNS is very CPU intensive, it is crucial that the timings are as good as 
possible. The current setup has been tested on a single processor with different 
chemistry and transport data, and the results are presented in Table [21 Run A, with 
the full Li mechanism and mixture averaged transport coefficients, use the most 
resources, as expected. By simplifying the transport data [21] (Run C) Eq. (flTI ) is 
substituted by 

(33) 



D, = Dn — 



and Eq. (|23l) is substituted by 

A=pCpKoT" (34) 

where n = 0.7 and Dq = kq = 2.9 x 10"^ g/(s cm K") leading to a 28% reduction in 
CPU consumption. Lets now turn off reactions, but still keeping all the 13 species 
(Run D), and an additional 53% reduction is achieved. 

For comparison. Run E is shown in order to see how much is gained by solving 
only the Navier-Stokes equation together with the continuity equation, assuming 
an isothermal medium with transport coefficients and thermodynamics such that 
all species can be neglected. It is seen that this is 20 times faster than Run D. This 
large difference is due to the fact that for Run E only 4 equations are solved, in 
contrast to the 18 equations for run D. Furthermore, and even more importantly, 
the time consuming process of determining the thermodynamics, such as enthalpy 
and heat capacity, together with the calculation of the viscosity, is omitted. 
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Figure 7: 3D images of the hot spot at moment t 
panel). 



s {left panel) and r=1.2xl0s {right 



5.2. Spherical flame front 

A study of the spherical and cylindrical flames is important because these 
cases are useful for determining important parameters in premixed combustion 
such as burning velocity, flame stretch rate, and flame curvature. There is a lot of 
numerical and experimental research in this area [|l3l] . The most important diffi- 
culty in the numerical approach is the large computational demand. The typical 
mesh size has to be 6x = AO- 60 /im [141 . i.e. for a cube of 3 cm"* one needs about 
500^ grid points and ideally 512 processors. 

For illustration purposes we consider a smaller cube (1 cm^), centered at the 
reference point with the hot spherical spot in its center (see Fig. [8]). The initial 
hydrogen-air mixture with Y^-, = 2.4 %, Yq^ = 23 % and Fn2 = 74.6 % is under a 
pressure of ;? = 1 bar. We use NSCBC boundary conditions, take a grid size of (80 
X 80 X 80), and 25 processors on the Cray XT4/XT5. The results are presented 
in Figs. |7] and [81 The 3D images of the hot spot at the diff"erent moments (at t = 
s and t = 1.2 X 10 "* s) are presented in Fig. Ul In Fig. [8] one sees that the gas 
is burned in the center and then the flame front is expanding symmetrically in all 
three directions. 

This problem is also used as a good test for the fully three dimensional NSCBC 
boundary conditions. We tested the implemented NSCBC boundary condition 
both for laminar and turbulent regimes. In the laminar case we find that due to 
the full NSCBC boundary conditions nlaA the code runs well up to the moment 
when the flame front comes to the domain boundaries. In the turbulent regime the 
problems appear near the corners and edges of the domain because of the eddies 
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Figure 8: Temperature as a function of x coordinate in the mid-plane of the box at f = s (dashed 
curve), t — 10^^ s {dotted-dashed curve) and t - 1.9 x 10"'* s {solid curve). 



at the boundaries. We avoid such a problem by using buffer (or sponge) zones (for 
details see tM)- We add the term to the right-hand side of the momentum equation 






V! - Vref,, 



^(Xi), j=l,...,Ni, 



(35) 



where j denotes the meshpoint and A^, is total number of grid points in the i di- 
rection, and dots indicate the presence of terms that where already specified in 
equation Q, ^(x,) is equal to zero everywhere except in the buffer zones where 
it is equal to unity. The length of the buffer zone is 10 % of the domain, and we 
choose T = 56t and Vief,; = 0. 

6. Conclusions 

In this paper we have presented a high-order public domain code for direct 
numerical simulation of compressible flows with detailed chemical reactions. The 
Pencil Code provides sixth-order spatial accuracy in the simple one-step reaction 
case, and fifth order accuracy in the case where upwinding for density advection is 
necessary. For validation purposes we compare our results with the Chemkin tool 
for OD and ID test problems, and show that they are in good agreement. Finally, 
we calculate the flame speed in 3D both in laminar and turbulent cases. 
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The code is well suited for considering also more complicated reaction schemes 
such as methane combustion. Furthermore, it is straightforward to consider the 
interaction with additional chemicals such as nitrogen and to follow the produc- 
tion of NOx gases. In particular, it is important to consider combustion in the 
presence of steam. This is well known to lead to a reduction of NOx gases. 
Combustion in the presence of more complicated boundary conditions involv- 
ing, for example, smaller inlet geometries has also been considered. Some of 
these cases, including those with a turbulent inlet, are available among the many 
sample cases that come with the code. For the benefit of the community, it is 
advantageous if prospective contributers to the code ask one of the code owners 
listed on http : //pencil-code . googlecode . c6m7| to obtain permission as a 
committer. 
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